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We study the pseudogaps in the spectra of the half-filled 2D 
Hubbard model using both finite-size and dynamical cluster 
approximation (DCA) quantum Monte Carlo calculations. A 
charge pseudogap, accompanied by non-Fermi liquid behavior 
in the self energy, is shown to persist in the thermodynamic 
limit. The DCA (finite-size) method systematically underes- 
timates (overestimates) the width of the pseudogap. A spin 
pseudogap is not seen at half-filling. 



Introduction For over a decade it has been recognized 
that the normal state properties of high-T c superconduc- 
tors are unusual and appear to have non-Fermi liquid 
characteristics. M One of the most remarkable features 
of the normal state is a suppression of the density of 
states at the Fermi energy in a temperature regime above 
T c in underdoped samples. Angular resolved photoemis- 
sion experiments [|2|,[3| show that this pseudogap in the 
spectral function has a d-wave anisotropy, the same sym- 
metry as the superconducting order parameter in these 
materials. This, along with theories that short-ranged 
spin fluctuations mediate pairing in the high-T c cuprates 
pLpl, emphasizes the importance of understanding the 
normal state, insulating phase. 

It is thought by many that the two-dimensional Hub- 
bard model, or closely related models, should capture the 
essential physics of the high-T c cuprates. M Yet, despite 
years of effort, neither the precursor pseudogap nor d- 
wave superconducting order have been conclusively seen 
in the Hubbard model. 

Intuitively, one may expect that the Hubbard model 
should show pseudogap behavior. At half-filling, the 
ground state of the 2D Hubbard model is an antifer- 
romagnetic insulator |qj7| and the spectrum is there- 
fore gapped. However, the Mermin- Wagner theorem pre- 
cludes any transition at finite T, so as the temperature 
is lowered one may anticipate that a pseudogap will de- 
velop. |§] This question has been previously addressed in 
the 2D Hubbard insulator by finite-size lattice Quantum 
Monte Carlo (QMC) |j 10 and approximate many-body 
techniques [[ll|-|l3|. The results have been contradictory 
and inconclusive as to the existence of a pseudogap at 
low temperatures, due to limitations of these techniques. 

Using the recently developed Dynamical Cluster Ap- 
proximation (DCA) 0,^5) we find that at sufficiently 
low temperatures a pseudogap opens in the single par- 
ticle spectral weight A(k, ui) of the 2D Hubbard model 



with a simultaneous destruction of the Fermi liquid state 
due to critical fluctuations above the T = transition 
temperature. This occurs in the weak-to-intermediate 
coupling regime U < W, where U is the on-site Coulomb 
energy and W the non-interacting band width. 

Using finite-sized techniques, it is difficult to determine 
if a gap persists in the thermodynamic limit. At half fill- 
ing, finite-size QMC calculations display a gap in their 
spectra as soon as the correlation length exceeds the lat- 
tice size, so they tend to overestimate the pseudogap as 
it would appear in the thermodynamic limit. Finite-size 
scaling is complicated by the lack of an exact scaling 
ansatz for the gap and the cost of performing simula- 
tions of large systems. Calculations employing Dynami- 
cal Mean Field Approximation (DMFA) [|16| in the para- 
magnetic phase do not display this behavior since they 
take place in the thermodynamic limit rather than on 
a finite-size lattice. However, the DMFA lacks the non- 
local spin fluctuations often believed to be responsible 
for the pseudogap. The Dynamical Cluster Approxima- 
tion (DCA) is a fully causal approach which systemati- 
cally incorporates non-local corrections to the DMFA by 
mapping the problem onto an embedded impurity clus- 
ter of size N c . N c determines the order of the approx- 
imation and provides a systematic expansion parameter 
1/N C . While the DCA becomes exact in the limit of large 
N c it reduces to the DMFA for N c = 1. Thus, the DCA 
differs from the usual finite size lattice calculations in 
that it is a reasonable approximation to the lattice prob- 
lem even for a "cluster" of a single site. Like the DMFA, 
the DCA solution remains in the thermodynamic limit, 
but the dynamical correlation length is restricted to the 
size of the embedded cluster. Thus the DCA tends to 
underestimate the pseudogap. 

Method The DCA is based on the assumption that 
the lattice self energy is weakly momentum dependent. 
This is equivalent to assuming that the dynamical inter- 
site correlations have a short spatial range b < L/2 where 
L is the linear dimension of the cluster. Then, according 
to Nyquist's sampling theorem jTj]], to reproduce these 
correlations in the self energy, we only need to sample the 
reciprocal space at intervals of Afc « 2ir/L. Therefore, 
we could approximate G(K + k) by G(K) within the cell 
of size (tt/L) (see, Fig. [j]) centered on the cluster mo- 
mentum K (wherever feasible, we suppress the frequency 
labels) and use this Green function to calculate the self 
energy. Knowledge of these Green functions on a finer 
scale in momentum is unnecessary, and may be discarded 
to reduce the complexity of the problem. Thus the clus- 



ter self energy can be constructed from the coarse-grained 
average of the single-particle Green function within the 
cell centered on the cluster momenta: 



N, 



G(K) = -^G(K + k), 



(1) 



where N is the number of points of the lattice, N c is 
the number of cells in the cluster, and the k summa- 
tion runs over the momenta of the cell about the clus- 
ter momentum K (see, Fig. [j]). For short distances 
r < L/2 the Fourier transform of the Green function 
G(r) w G(r) + 0((YAfc) 2 ), so that short ranged corre- 
lations are reflected in the irreducible quantities con- 
structed from G; whereas, longer ranged correlations 
r > L/2 are cut off by the finite size of the cluster. [[Uj 




FIG. 1. N c = 4 cluster cells (shown by different fill pat- 
terns) that partition the first Brillouin Zone (dashed line). 
Each cell is centered on a cluster momentum K (filled circles) . 
To construct the DCA cluster, we map a generic momentum 
in the zone such as k to the nearest cluster point K = M(k) 
so that k = k — K remains in the cell around K. 

This coarse graining procedure and the relationship of 
the DCA to the DMFA is illustrated by a microscopic 
diagrammatic derivation of the DCA. For Hubbard-likc 
models, the properties of the bare vertex are completely 
characterized by the Laue function A which expresses 
the momentum conservation at each vertex. In a con- 
ventional diagrammatic approach A(ki, k 2 ,k 3 , k 4 ) = 
J2 r ex P [* r ' ( k i - k 2 + k 3 - k 4 )] = iV<5 kl+k2!k3+k4 where 
ki and k 2 (k 3 and k 4 ) are the momenta entering (leaving) 
each vertex through its legs of G. However as D — > oo 
Miillcr-Hartmann showed that the Laue function reduces 
to [IH 



A D ^ oc (k 1 , k 2 , k 3 , k 4 ) = 1 + 0(1/ D) 



(2) 



The DMFA assumes the 

same Laue function, Admf A (ki > ^ , k 3 , k 4 ) = 1, even in 



the context of finite dimensions. Thus, the conservation 
of momentum at internal vertices is neglected. Therefore 
we may freely sum over the internal momenta at each 
vertex in the generating functional &dmfa- This leads 
to a collapse of the momentum dependent contributions 
to the functional &dmfa and only local terms remain. 

The DCA systematically restores the momentum con- 
servation at internal vertices. As discussed above the 
Brillouin-zone is divided into N c = L D cells of size 
(2n/L) D . Each cell is represented by a cluster momen- 
tum K in the center of the cell. We require that mo- 
mentum conservation is (partially) observed for momen- 
tum transfers between cells, i.e. for momentum transfers 
larger than Afc = 2?t/L, but neglected for momentum 
transfers within a cell, i.e less than Afc. This require- 
ment can be established by using the Laue function |14J| 



A_DCA(ki,k 2 ,k 3 ,k 4 



N r 5 



M(ki)+M(k 3 ),M(k 2 )+M(k 4 ) 



(3) 



where M(k) is a function which maps k onto the mo- 
mentum label K of the cell containing k (see, Fig. [y). 

With this choice of the Laue function the momenta 
of each internal leg may be freely summed over the cell. 
This is illustrated for the second-order term in the gener- 
ating functional in Fig. @. Thus, each internal leg G(ki) 
in a diagram is replaced by the coarse-grained Green 
function G(M(ki)), defined by Eq. | 
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M(k,) +M(k 2 ) ,M(k,) +M(k 4 ) 



f EG(K+q) = 



G(K) 




K'+Q 



FIG. 2. The DCA choice of the Laue function Eq. H leads 
to the replacement of the lattice propagators G(ki), G(ka),... 
by coarse grained propagators G(K), G(K'), ... (Eq. hi) in the 
internal legs of <&dca, illustrated for a second order diagram. 

The diagrammatic sequences for the generating func- 
tional and its derivatives are unchanged; however, the 
complexity of the problem is greatly reduced since iV c <c 
N. We showed previously @,||| that the DCA estimate 
of the lattice free-energy is minimized by the approxima- 
tion E(k) « S(M(k)), where 5<P DCA /5G = S. 

The cluster problem is then solved by usual techniques 
such as QMC [ ll| , the non-crossing approximation |lj| or 
the Fluctuation-Exchange approximation. Here we em- 
ploy a generalization of the Hirsh-Fye QMC algorithm 
pOf to solve the cluster problem. The initial Green func- 
tion for this procedure is the bare cluster Green function 
C(K)- 1 = G(K)^ 1 + £(K) which must be introduced to 
avoid over-counting diagrams. The QMC estimate of the 
cluster self energy is then used to calculate a new esti- 
mate of G(K) using Eq. nl The corresponding G(K) is 
used to reinitialize the procedure which continues until 
the self energy converges to the desired accuracy. 




FIG. 3. The spectral density A(k, u), and the real 
Re£(k, ui) and imaginary ImE(k, to) parts of the self-energy 
for the 2D Hubbard model via the DCA with a paramagnetic 
host at k = (-7T, 0) for a 64-site cluster (N c = 64) at various 
temperatures. The on-site Coulomb repulsion U — 5.2, the 
band width W — 8, and the filling (n) = 1. As the tem- 
perature is lowered, the system first builds a Fermi-liquid-like 
peak in A(k,ui). By (5 — 2.6, a pseudogap begins to develop 
in A(k,u>) and simultaneously, ReE(k,o;) develops a positive 
slope at to — 0, a signal of a non-Fermi liquid. The pseudogap 
deepens as the temperature is further lowered. The imaginary 
part of the dynamic spin susceptibility, divided by the static 
spin susceptibility is shown in the inset. No spin gap is seen. 

Results We study the 2D Hubbard Hamiltonian: 

+U^2{n^ - -){n ti - -) -fiJ2 n ^- ( 4 ) 



where c\ a (cia) creates (destroys) an electron at site i with 
spin a, U is the on-site Coulomb potential, and rii a = 
c icr c i<? i s ^ ne number operator. We set the overlap integral 
t = 1 and measure all energies in terms of t. We work 
at /i = where the system is half- filled ((n) — I). We 
choose U = 5.2, which is well below the value U > W 
believed to be necessary to open a Mott-Hubbard gap. 

We also calculate the angle integrated dynamical spin 
susceptibility shown in the inset. It does not have a 
pseudogap, as expected for the half-filled model since the 
spin-wave spectrum is gapless. Since a (spin) charge gap 
is generally defined as one which appears in the (spin) 
charge dynamics or thermodynamics, we conclude that 
the pseudogap is only in the charge response and is due 
to short-ranged antiferromagnetic spin correlations. 

Fig. H shows the spectral density A(k, u>), and the 
real Re£(k, uj) and imaginary Im£(k, uj) parts of the 
self-energy for the 2D Hubbard model via the DCA 
with a paramagnetic host at the Fermi surface X point 



k = (it, 0) for a 64-site cluster (N c — 64) at various 
temperatures. We obtain the spectral function A(k, uj) 
via the Maximum Entropy Method (MEM). Q As the 
temperature is lowered, the system first builds a Fermi- 
liquid-like peak in A(k, to). By j3 = 2.6, a pseudogap 
begins to develop in A(k, ui). The pseudogap builds as 
the temperature is further lowered. 

Fig. shows the spectral function A(k, uj) at the half- 
filled Fermi surface point k = (tt/2, 7t/2). The qualitative 
features are similar to Fig. S, but the pseudogap opens at 
a lower temperature and the distance between the peaks 
is less than that seen at the X point. This behavior is 
reminiscent of the anisotropy of the pseudogap observed 
experimentally in the insulating [|3| and in the supercon- 
ducting g cuprates, but is not large enough to be compa- 
rable with that seen experimentally. We speculate that 
the anisotropy seen here may be due to a difference in 
the number of states near the Fermi energy at these two 
points in the zone. 




FIG. 4. The spectral density A(k,ui), and the real 
ReE(k, ui) and imaginary ImE(k, to) parts of the self-energy 
for the 2D Hubbard model via the DCA at k = (tt/2, tt/2) 
and the same parameters as Fig. H. Again the system first 
builds a Fermi-liquid-like peak in A(k, uj) and then develops 
a pseudogap in A(k,uj) with a simultaneous non-Fermi liq- 
uid behavior in ReE(k,o;). Here, though, the pseudogap first 
appears at a lower temperature than at k = (n,0). 

The DCA self-energy spectra Figs. S & support the 
spectral evidence. At the X point, the slope of the real 
part Re£(k, u>) becomes positive below (3 = 2.6, the tem- 
perature at which we observed the opening of a pseudo- 
gap. This signals the appearance of two new solutions in 
the quasiparticle equation Re(w — ek — S(k, uj)) = in 
addition to the strongly damped solution at uj = which 
is also present in the noninteracting system. These two 
new quasiparticle solutions for the same k-vector indi- 
cate precursor effects of the onset of antiferromagnetic 
ordering which entails a doubling of the unit cell. They 



are referred as shadow states and are caused by antiferro- 
magnetic spin fluctuations in the paramagnetic state. At 
these temperatures, the imaginary part Im£(k,w) dis- 
plays a local minimum at u> = indicating the break- 
down of the Fermi liquid behavior. We note that a differ- 
ent conclusion was previously reached in a FLEX study 
jflj, which found that ImS(k, u>) has a local minimum 
at u) = which was not accompanied by an opening of 
a pseudogap. Since the pseudogap is due to short-range 
spin correlations, we conclude that FLEX underestimates 
these correlations. 
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FIG. 5. The spectral density A(k,cj) at k = (tt, 0) for the 
2D Hubbard model via (a) the DCA and (b) finite-size Quan- 
tum Monte Carlo (QMC) at an inverse temperature times the 
bandwidth fiW = 40 on various size clusters. The tempera- 
ture T* at which the pseudogap first becomes apparent in the 
DCA spectra, as well as the full width A measured from peak 
to peak is plotted in the inset. The finite-size QMC over- 
estimates A and T* , whereas the DCA QMC systematically 
underestimates them. 

It is instructive to compare the DCA results with those 
obtained by finite-size QMC calculations. Fig. g shows 
the spectral density A(tv,0,u>) obtained by analytically 
continuing both finite-size and DCA QMC data. In spite 
of the difference in the two methods, the information they 
provide is complimentary. In the finite-size results, (b), 
we see a similar opening of a pseudogap. However, as the 
length of the antifcrromagnetic (AF) correlations reach 
the longest length on the finite-size lattice, the system 
develops a full gap. Thus, the finite-size QMC overesti- 
mates the size of the gap. In the DCA results, (a), the 
pseudogap emerges as soon as N c > 1. The temperature 
T* at which the pseudogap first becomes apparent in the 
spectra, as well as the full width A measured from peak 
to peak is plotted in the inset. Both T* and A increase 
with N c . Since the DCA calculation remains in the ther- 
modynamic limit, a full gap due to antiferromagnetic cor- 
relations alone cannot open until their correlation length 
diverges. However, since these correlations are restricted 



to the size of the cluster, the DCA systematically under- 
estimates the size of the gap. Thus, if a pseudogap exists 
in the DCA for finite N c , it should persist in the limit as 
N c -> oo. 

In summary, we have employed the recently developed 
DCA to study the long-open question of whether the half- 
filled Hubbard model has a pseudogap due to AF spin 
fluctuations. We find conclusive evidence of a pseudogap 
in the charge dynamics and have shown unambiguously 
that the T = phase transition of the half-filled model 
is preceded by an opening of a pseudogap in A(k.F,cj) 
accompanied by pronounced non-Fermi liquid behavior 
in E(kp>, uj). 
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